可以使用Python的Sympy库求微分方程的符号解，Scipy库求微分方程的数值解;本文介绍数值解。
数值解¶
odeint¶
Python 对常微分方程的数值求解是基于一阶方程进行的, 高阶微分方程必须化成一阶方程组, 通常采用龙格—库塔方法。scipy.integrate 模块的 odeint 函数求常微分方程的数值解, 其基本调用格式为:
$$
\text { sol=odeint(func, y0, t) }
$$
其中 func 是定义微分方程的函数或匿名函数, y0 是初始条件的序列, $\mathbf{t}$ 是 一个自变量取值的序列 (t 的第一个元素一定为初始时刻), 返回值 sol 是 对应于序列 $\mathrm{t}$ 中元素的数值解, 如果微分方程组中有 $n$ 个函数, 返回值 sol 是 $n$ 列的矩阵, 第 $i(i=1,2, \cdots, n)$ 列对应于第 $i$ 个函数的数值解。
例1¶
求微分方程 $$ \left\{\begin{array}{l} y^{\prime}=-2 y+x^{2}+2 x, \\ y(1)=2 . \end{array}\right. $$ 在 $1 \leq x \leq 10$ 步长间隔为 $0.5$ 点上的数值解。
代码
from scipy.integrate import odeint
from numpy import arange
dy=lambda y, x: -2*y+x**2+2*x
x=arange(1, 10.5, 0.5)
sol=odeint(dy, 2, x)
print("x={}\n对应的数值解y={}".format(x, sol.T))
例2¶
求下述微分方程的解 $$ \left\{\begin{array}{l} \frac{d^{2} y}{d x^{2}}+2 \frac{d y}{d x}+2 y=0 \\ y(0)=0, \quad y^{\prime}(0)=1 \end{array}\right. $$
解: 引进 $y_{1}=y, y_{2}=y^{\prime}$, 则可以把原来的二阶微分方程化为如下一阶 微分方程组: $$ \begin{cases}y_{1}^{\prime}=y_{2}, & y_{1}(0)=0, \\ y_{2}^{\prime}=-2 y_{1}-2 y_{2}, & y_{2}(0)=1 .\end{cases} $$ 求数值解和画图的程序如下:
from scipy.integrate import odeint
from sympy.abc import t
import numpy as np
import matplotlib.pyplot as plt
def Pfun(y,x):
y1, y2=y;
return np.array([y2, -2*y1-2*y2])
x=np.arange(0, 10, 0.1) #创建时间点
sol1=odeint(Pfun, [0.0, 1.0], x) #求数值解
plt.rc('font',size=16); plt.rc('font',family='SimHei')
plt.plot(x, sol1[:,0],'r*',label="数值解")
plt.plot(x, np.exp(-x)*np.sin(x), 'g', label="符号解曲线")
plt.legend()
plt.show()
例3¶
Lorenz 模型的混沌效应。
Lorenz 模型是由美国气象学家 Lorenz 在研究大气运动时, 通过简化对流模型, 只保留 3 个变量提出的一个完全确定性的一阶自治常微分方程组 (不显含时间变量), 其方程为 $$ \left\{\begin{array}{l} \dot{x}=\sigma(y-x), \\ \dot{y}=\rho x-y-x z, \\ \dot{z}=x y-\beta z . \end{array}\right. $$ 其中, 参数 $\sigma$ 为 Prandtl 数, $\rho$ 为 Rayleigh 数, $\beta$ 为方向比。
代码
from scipy.integrate import odeint
import numpy as np
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
def lorenz(w,t):
sigma=10; rho=28; beta=8/3
x, y, z=w;
return np.array([sigma*(y-x), rho*x-y-x*z, x*y-beta*z])
t=np.arange(0, 50, 0.01) #创建时间点
sol1=odeint(lorenz, [0.0, 1.0, 0.0], t) #第一个初值问题求解
sol2=odeint(lorenz, [0.0, 1.0001, 0.0], t) #第二个初值问题求解
plt.rc('font',size=16); plt.rc('text',usetex=True)
ax1=plt.subplot(121,projection='3d')
ax1.plot(sol1[:,0], sol1[:,1], sol1[:,2],'r')
ax1.set_xlabel(r'$x$'); ax1.set_ylabel(r'$y$'); ax1.set_zlabel(r'$z$')
ax2=plt.subplot(122,projection='3d')
ax2.plot(sol1[:,0]-sol2[:,0], sol1[:,1]-sol2[:,1], sol1[:,2]-sol2[:,2],'g')
ax2.set_xlabel(r'$x$'); ax2.set_ylabel(r'$y$'); ax2.set_zlabel(r'$z$')
plt.tight_layout()
plt.savefig("images/diff0602.png", dpi=500);
plt.show()
print("sol1=",sol1, '\n\n', "sol1-sol2=", sol1-sol2)
Lorenz 模型如今已经成为混沌领域的经典模型, 第一个混沌吸引子一 Lorenz 吸引子也是在这个系统中被发现的。系统中三个参数的选择对系统 会不会进入混沌状态起着重要的作用。下图 给出了 Lorenz 模型在 $\sigma=10, \rho=28, \beta=8 / 3$ 时系统的三维演化轨迹。由图 可见, 经 过长时间运行后, 系统只在三维空间的一个有限区域内运动, 即在三维相空 间里的测度为零。图 (A) 显示出我们经常听到的“蝴蝶效应”。图 (B) 给出了系统从两个靠的很近的初值出发（相差仅 0.0001）后, 解的偏差演化曲线。
随着时间的增大，可以看到两个解的差异越来越大，这正是动力学系统对初值敏感性的直观表现，由此可断定此系统的这种状态为混沌态。混沌运动是确定性系统中存在随机性，它的运动轨道对初始条件极端敏感。
